

/** @file basic_sector_decomposition.cc
 *
 * Implementation of basic_sector_decomposition 
 *
 */

/*
 *  Copyright (C) 2007 Stefan Weinzierl
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include <vector>
#include <iostream>
#include <ctime>
#include <sstream>


#include "basic_sector_decomposition.h"
#include "monte_carlo.h"
#include "constants.h"

namespace sector_decomposition {

  /**
   *
   * The main routine to compute an integral of the form 
   * \f[
   *     \int\limits_{x_j \ge 0} d^nx \;\delta(1-\sum_{i=1}^n x_i)
   *      \left( \prod\limits_{i=1}^n x_i^{a_i+\varepsilon b_i} \right)
   *       \prod\limits_{j=1}^r \left[ P_j(x) \right]^{c_j+\varepsilon d_j}.
   * \f]
   * The integration is over the standard simplex:
   * \f[
         x_j \ge 0, \;\;\; x_1+...+x_n=1.
   * \f]
   * For integration over the standard hyper-cube, the function do_sector_decomposition_cube can be used.
   * 
   * The arguments are as follows:
   * \verbatim
       monte_carlo_result do_sector_decomposition(const integration_data & global_data,
                                                  const integrand & integrand_in,
                                                  const monte_carlo_parameters & mc_parameters,
                                                  int verbose_level = 0);
     \endverbatim
   * The input are three structures, integration_data, integrand,
   * monte_carlo_parameters, which are described in detail below, and an optional parameter verbose_level.
   * With the help of the optional parameter verbose_level one can choose the amount 
   * of information the program prints out during the run.
   * The function returns a structure monte_carlo_result, which again is described below.
   *
   * The class integration_data contains the data which will not be modified
   * by the algorithm for the sector decomposition.
   * It has a constructor of the form
   * \verbatim
      integration_data(const std::vector<GiNaC::ex> & list_feynman_parameter, 
                       GiNaC::ex epsilon, int order);
     \endverbatim
   * where list_feynman_parameter is a vector holding the symbols of the Feynman parameters,
   * epsilon is the symbol of the expansion parameter
   * and order defines which term of the Laurent series in \f$\varepsilon\f$ should be computed.
   * 
   * The integrand
   * \f[
   *     \left( \prod\limits_{i=1}^n x_i^{a_i+\varepsilon b_i} \right)
   *      \prod\limits_{j=1}^r \left[ P_j(x) \right]^{c_j+\varepsilon d_j}
   * \f]
   * is encoded in the class integrand. This class has a constructor of the form
   * \verbatim
       integrand(const std::vector<exponent> & nu,
                 const std::vector<GiNaC::ex> & poly_list, 
                 const std::vector<exponent> & c);
     \endverbatim
   * where nu is a vector of size \f$n\f$ holding the exponents \f$a_i+\varepsilon b_i\f$.
   * poly_list is a vector of size \f$r\f$, holding the polynomials \f$P_j\f$ in the Feynman parameters.
   * The corresponding exponents are stored in the vector c, again of size \f$r\f$.
   * As exponents are generally of the form \f$a+b\varepsilon\f$, a special class is available for them with the constructor
   * \verbatim
      exponent(int a, int b);
     \endverbatim
   * Here \f$a\f$ and \f$b\f$ are integer numbers.
   * Rational numbers for \f$a\f$ and \f$b\f$ are also allowed, then the constructor
   * \verbatim
      exponent(GiNaC::ex order_0, GiNaC::ex order_1);
     \endverbatim
   * has to be used. In this case order_0 and order_1 have to be of type GiNaC::numeric. 
   *
   * In applications one encounters often integrands where the explicit powers of the 
   * Feynman parameters are missing. For integrands of the form
   * \f[
        \prod\limits_{j=1}^r \left[ P_j(x) \right]^{c_j+\varepsilon d_j}
   * \f]
   * there is a simpler constructor of the form
   * \verbatim
      integrand(size_t n, const std::vector<GiNaC::ex> & poly_list, 
                const std::vector<exponent> & c);
     \endverbatim
   * Here, \f$n\f$ is the number of Feynman parameters.
   *
   * Parameters associated to the Monte Carlo integration are specified with the help of the class
   * monte_carlo_parameters.
   * This class is constructed as follows:
   * \verbatim
      monte_carlo_parameters(size_t iterations_low, size_t iterations_high, 
                             size_t calls_low, size_t calls_high);
     \endverbatim
   * The program uses the Vegas-algorithm based on an adaptive grid. 
   * The program does first a Monte Carlo integration with iterations_low iterations
   * with calls_low function evaluations each. This phase is solely used to adapt the grid.
   * The numerical result of the Monte Carlo integration is then obtained from the
   * second stage with iterations_high iterations of calls_high function calls each.
   *
   * The main function do_sector_decomposition returns the numerical results
   * of the Monte Carlo integration in the class monte_carlo_result.
   * This class has the methods
   * \verbatim
      class monte_carlo_result {
       public :
        double get_mean(void) const;
        double get_error(void) const;
        double get_chi_squared(void) const;
      };
     \endverbatim
   * which return the mean value of the Monte Carlo integration, the error estimate and the associated \f$\chi^2\f$.
   *
   * The strategy which is used for the iterated sector decomposition can be selected with the global variable
   * sector_decomposition::CHOICE_STRATEGY. Available keywords are STRATEGY_A, STRATEGY_B, STRATEGY_C, STRATEGY_X
   * and STRATEGY_Y. The strategies A, B and C are guaranteed to terminate.
   * Strategies X and Y are heuristic strategies, for which there is no proof that they terminate.
   * They are included, since in the cases where they terminate the number of generated sub-sectors tends to be smaller
   * than the corresponding numbers for the other strategies.
   * The default is strategy C. 
   *
   * The function do_sector_decomposition is just the entry point and calls the function do_sector_decomposition_simplex, which
   * in turn delegates evaluation to  a function object
   * of type basic_sector_decomposition.
   *
   */
  monte_carlo_result do_sector_decomposition(const integration_data & global_data,
					     const integrand & integrand_in, 
					     const monte_carlo_parameters & mc_parameters,
					     int verbose_level)
  {
    return do_sector_decomposition_simplex(global_data, integrand_in, mc_parameters, verbose_level);
  }

  /**
   *
   * Entry point for the integration over the standard simplex.
   *
   * The function do_sector_decomposition calls do_sector_decomposition_simplex.
   *
   * The arguments of the function are as for the function do_sector_decomposition.
   *
   */
  monte_carlo_result do_sector_decomposition_simplex(const integration_data & global_data,
						     const integrand & integrand_in, 
						     const monte_carlo_parameters & mc_parameters,
						     int verbose_level)
  {
    basic_sector_decomposition T;
    return T(global_data,integrand_in,mc_parameters,SIMPLEX,-1);
  }

  /**
   *
   * Entry point for the integration over the standard hyper-cube.
   *
   * This routine computes an integral of the form 
   * \f[
   *     \int\limits_0^1 d^nx \;
   *      \left( \prod\limits_{i=1}^n x_i^{a_i+\varepsilon b_i} \right)
   *       \prod\limits_{j=1}^r \left[ P_j(x) \right]^{c_j+\varepsilon d_j}.
   * \f]
   * The integration is over the standard hyper-cube:
   * \f[
         0 \le x_j \le 1.
   * \f]
   * The polynomials \f$P_j(x)\f$ are assumed to be non-vanishing inside the hyper-cube,
   * but may vanish on the boundary of the hyper-cube.
   * This means that the polynomials may vanish, if one or more variables \f$x_j\f$ are equal
   * to zero or one.
   * 
   * If the polynomials vanish only, if one or more variables are zero, the funtion
   * do_sector_decomposition_cube_coordinate_subspace
   * can be used, which is more efficient for this case.
   *
   * The arguments of the function are as for the function do_sector_decomposition.
   *
   */
  monte_carlo_result do_sector_decomposition_cube(const integration_data & global_data,
						  const integrand & integrand_in, 
						  const monte_carlo_parameters & mc_parameters,
						  int verbose_level)
  {
    basic_sector_decomposition T;
    return T(global_data,integrand_in,mc_parameters,CUBE,verbose_level);
  }

  /**
   *
   * Entry point for the integration over the standard hyper-cube, with singularities only on coordinate subspaces.
   *
   * This routine computes an integral of the form 
   * \f[
   *     \int\limits_0^1 d^nx \;
   *      \left( \prod\limits_{i=1}^n x_i^{a_i+\varepsilon b_i} \right)
   *       \prod\limits_{j=1}^r \left[ P_j(x) \right]^{c_j+\varepsilon d_j}.
   * \f]
   * The integration is over the standard hyper-cube:
   * \f[
         0 \le x_j \le 1.
   * \f]
   * The polynomials \f$P_j(x)\f$ are assumed to be non-vanishing inside the hyper-cube,
   * but may vanish on the coordinate subspaces of the hyper-cube.
   * This means that the polynomials may vanish, if one or more variables \f$x_j\f$ are equal
   * to zero.
   * 
   * The arguments of the function are as for the function do_sector_decomposition.
   *
   */
  monte_carlo_result do_sector_decomposition_cube_coordinate_subspace(const integration_data & global_data,
								      const integrand & integrand_in, 
								      const monte_carlo_parameters & mc_parameters,
								      int verbose_level)
  {
    basic_sector_decomposition T;
    return T(global_data,integrand_in,mc_parameters,CUBE_COORDINATE_SUBSPACE,verbose_level);
  }

  /**
   *
   * The old entry point, which is kept for backwards compatibility.
   *
   */
  monte_carlo_result basic_sector_decomposition::operator() (const integration_data & global_data,
							     const integrand & integrand_in, 
							     const monte_carlo_parameters & mc_parameters,
							     const int & verbose_level)
  {
    return this->operator()(global_data, integrand_in, mc_parameters, SIMPLEX, verbose_level);
  }

  /**
   *
   * The main routine of the class basic_sector_decomposition.
   *
   * The routine performs the following steps:
   *
   * 1. If the integration domain is the standard simplex, convert all polynomials to homogeneous polynomials.
   *
   * 2. Decompose the integral into primary sectors.
   *
   * 3. Iteration of the sector decomposition.
   *
   * 4. Generate the subtraction terms.
   *
   * 5. Expand in \f$\varepsilon\f$.
   *
   * 6. Numerical integration.
   *
   */
  monte_carlo_result basic_sector_decomposition::operator() (const integration_data & global_data,
							     const integrand & integrand_in, 
							     const monte_carlo_parameters & mc_parameters,
							     const int & integration_domain,
							     const int & verbose_level)
  {
    time_t time_start;
    time_t time_end;

    if ( verbose_level>1 )
      {
	std::cout << "Input integrand:" << std::endl;
	std::cout << integrand_in.get_integrand(global_data) << std::endl;
	std::cout << std::endl;
      }

    // --------------------------------------------------------------
    // test for homogeneous polynomials
    time(&time_start);
    if ( verbose_level>0 ) std::cout << "Homogeneous polynomials:" << std::flush;

    integrand integrand_homogeneous;

    if ( integration_domain == SIMPLEX )
      {
	generate_homogenous_polynomials(global_data, integrand_in, integrand_homogeneous);
      }
    else // CUBE, CUBE_COORDINATE_SUBSPACE
      {
	// no homogeneous polynomials needed
	integrand_homogeneous = integrand_in;
      }

    time(&time_end);
    if ( verbose_level>0 ) std::cout << " Done. Time in seconds: " << difftime(time_end,time_start) << std::endl;

    if ( verbose_level>1 )
      {
	std::cout << integrand_homogeneous.get_integrand(global_data) << std::endl;
	std::cout << std::endl;
      }

    // --------------------------------------------------------------
    // primary sectors
    time(&time_start);
    if ( verbose_level>0 ) std::cout << "Primary sectors:        " << std::flush;

    std::vector<integrand> primary_sectors;

    if ( integration_domain == SIMPLEX )
      {
	generate_primary_sectors(global_data, integrand_homogeneous, primary_sectors);
      }
    else if ( integration_domain == CUBE )
      {
	generate_primary_sectors_cube(global_data, integrand_homogeneous, primary_sectors);
      }
    else if ( integration_domain == CUBE_COORDINATE_SUBSPACE )
      {
	// no generation of primary sectors needed
	primary_sectors.push_back( integrand_homogeneous );
      }

    time(&time_end);
    if ( verbose_level>0 ) std::cout << " Done. Time in seconds: " << difftime(time_end,time_start) << std::endl;
    if ( verbose_level>0 ) std::cout << " " << primary_sectors.size() << " primary sectors." << std::endl;

    if ( verbose_level>1 )
      {
	for (int k=0; k<primary_sectors.size(); k++)
	  {
	    //std::cout << k << ": " << (primary_sectors[k]).get_integrand(global_data) << std::endl;
	  }
	std::cout << std::endl;
      }

    // --------------------------------------------------------------
    // iteration
    time(&time_start);
    if ( verbose_level>0 ) std::cout << "Iterated sectors:       " << std::flush;

    std::vector<integrand> sectors_out;

    iterate_sector_decomposition(global_data, primary_sectors, sectors_out);
                        std::cout<< std::endl;

    std::cout << "\n\nvalues_strings.append(\"\"\"";
    for (int kk=0; kk<sectors_out.size(); kk++) {

        std::vector<integrand> sectors_out2;
        sectors_out2.push_back(sectors_out[kk]);

        std::vector<GiNaC::ex> subtracted_terms;

        generate_subtraction_terms(global_data, sectors_out2, subtracted_terms);

        GiNaC::ex expansion_out;


        expand_in_epsilon(global_data, subtracted_terms, expansion_out);

            monte_carlo_result res1;
            numerical_integration(global_data, mc_parameters,expansion_out, res1, verbose_level);
            std::cout <<  res1.get_mean() << ", " << res1.get_error() << std::endl;
    }
    std::cout << "\"\"\")";


    std::vector<GiNaC::ex> subtracted_terms;

    generate_subtraction_terms(global_data, sectors_out, subtracted_terms);

    GiNaC::ex expansion_out;

    expand_in_epsilon(global_data, subtracted_terms, expansion_out);

    monte_carlo_result res;
    numerical_integration(global_data, mc_parameters, expansion_out, res, verbose_level);


    if ( verbose_level>0 )
      {
	std::cout << "Order " << pow(global_data.epsilon,global_data.order) << ": " 
		  << res.get_mean() << " +/- " << res.get_error() << std::endl;
      }

    // --------------------------------------------------------------

    return res;
  }

  /**
   *
   * Ensures that all polynomials are homogeneous.
   *
   */
  void basic_sector_decomposition::generate_homogenous_polynomials(const integration_data & global_data, 
								   const integrand & integrand_in, integrand & integrand_out)
  {
    integrand_out = integrand_in;
    integrand_out.make_homogeneous(global_data);
  }

  /**
   *
   * Generates the primary sectors.
   *
   */
  void basic_sector_decomposition::generate_primary_sectors(const integration_data & global_data, 
							    const integrand & integrand_in, std::vector<integrand> & lst_integrand_out)
  {
    int n = global_data.n;

    for (int k=0; k<n; k++)
      {
	integrand sector_k = integrand_in;
	sector_k.restrict_to_primary_sector(k, global_data);
	lst_integrand_out.push_back( sector_k );
      }
  }

  /**
   *
   * Generates the primary sectors for the integration over the standard hyper-cube.
   *
   */
  void basic_sector_decomposition::generate_primary_sectors_cube(const integration_data & global_data, 
								 const integrand & integrand_in, std::vector<integrand> & lst_integrand_out)
  {
    int n = global_data.n;

    multi_index_counter index_set(2,n);
    for( index_set.init(); !index_set.overflow(); index_set++)
      {
	integrand sector_k = integrand_in;
	sector_k.restrict_to_primary_sector_cube(index_set, global_data);
	lst_integrand_out.push_back( sector_k );
      }

  }

  /**
   *
   * Iterated sector decomposition. The strategy is chosen according to the 
   * global variable CHOICE_STRATEGY.
   *
   */
  void basic_sector_decomposition::iterate_sector_decomposition(const integration_data & global_data, 
								std::vector<integrand> & lst_integrand_in, std::vector<integrand> & lst_integrand_out)
  {
    if ( CHOICE_STRATEGY == STRATEGY_A )
      {
	strategy_A(global_data, lst_integrand_in, lst_integrand_out);
	return;
      }
    else if ( CHOICE_STRATEGY == STRATEGY_B )
      {
	strategy_B(global_data, lst_integrand_in, lst_integrand_out);
	return;
      }
    else if ( CHOICE_STRATEGY == STRATEGY_C )
      {
	strategy_C(global_data, lst_integrand_in, lst_integrand_out);
	return;
      }
    else if ( CHOICE_STRATEGY == STRATEGY_X )
      {
	strategy_X(global_data, lst_integrand_in, lst_integrand_out);
	return;
      }
    else if ( CHOICE_STRATEGY == STRATEGY_Y )
      {
	strategy_Y(global_data, lst_integrand_in, lst_integrand_out);
	return;
      }
  }


  std::string getIndex(GiNaC::ex ex) {
    std::ostringstream s;
    s << ex;
    std::string str = s.str();
    return str.substr(1, str.length());
  }
  /**
   *
   * Iterated sector decomposition: Strategy A, Zeillinger's strategy.
   *
   */
  void basic_sector_decomposition::strategy_A(const integration_data & global_data, 
					      std::vector<integrand> & lst_integrand_in, std::vector<integrand> & lst_integrand_out)
  {

      std::vector<std::string> lines;
             for (int i = 0; i < lst_integrand_in.size(); i++) {
                  std::ostringstream s;
                  s << "\n(" << getIndex(global_data.list_feynman_parameter[i]);
                  s << ",(";
                  for (int j = 0; j < lst_integrand_in.size(); j++) {
                      if (i != j) {
                          s << getIndex(global_data.list_feynman_parameter[j]);
                           if (!(j == lst_integrand_in.size() - 2 && i == lst_integrand_in.size() - 1 )) {
                           if (!(j == lst_integrand_in.size() - 1)) {
                                                    s << ",";
                                                    }
                                                }
                      }

                  }
                  s << ",),)";
                  lines.push_back(s.str());
             }

    std::cout << "sector_strings.append(\"\"\"";

    while ( lst_integrand_in.size() > 0 ) 
      {
	bool flag_continue = true;




	// get the last integral
	integrand local_integrand = lst_integrand_in.back();
	lst_integrand_in.pop_back();

	std::string line = lines.back();
	lines.pop_back();

	size_t n = local_integrand.index.size();
	size_t r = local_integrand.polynomial_list.size();

	// ----------------------------------------------------
	// Zeillinger's strategy
	std::vector<size_t> index_set;

	if ( local_integrand.find_min_max(index_set,global_data) )
	  {
	    // create two new subsectors
	    for (int k=0; k<2; k++)
	      {
		integrand sub_integrand = local_integrand;

		sub_integrand.create_subsector(k, index_set, global_data);

        //std::cout << sub_integrand << std::endl;

        std::ostringstream os;
        int kk = k == 0 ? 1 : 0;
        os << line << ", (" << getIndex(global_data.list_feynman_parameter[sub_integrand.index[index_set[k]]]) <<
         ",(" << getIndex(global_data.list_feynman_parameter[sub_integrand.index[index_set[kk]]]) << ",),)";
        lines.push_back(os.str());
		lst_integrand_in.push_back(sub_integrand);
	      } // k=0; k<2;

	    flag_continue = false;
	  }

	if ( flag_continue )
	  {
	    // one variable factorizes
	    multi_index_ordered index_k(n,1);
	    for( index_k.init(); !index_k.overflow(); index_k++)
	      {
		for (int j=0; j<r; j++) 
		  {
		    if ( local_integrand.check_zero(j, index_k, global_data) ) 
		      {
			// can factor the k'th variable from the j'the polynomial
			size_t k = index_k[0];
			local_integrand.polynomial_list[j] = normal(local_integrand.polynomial_list[j]
								    /global_data.list_feynman_parameter[local_integrand.index[k]]);

			(local_integrand.weight[k]).sum_up( local_integrand.poly_exponent[j] );

			flag_continue = false;
			//std::cout << local_integrand;
		      } // check_zero
		  } // j=0; j<r;
	      } // loop over index_k

	    // some variables were factorized
	    if ( !flag_continue) {
	    lines.push_back(line);
	    lst_integrand_in.push_back(local_integrand);
	    } else {
	                //std::cout << line;
	    }

	  }

	// ----------------------------------------------------
	// no zeros

	if ( flag_continue ) {
	 lst_integrand_out.push_back(local_integrand);
	 std::cout << line;
	}
      }

         std::cout << "\n\"\"\")";
  }

  /**
   *
   * Iterated sector decomposition: Strategy B, Spivakovsky's strategy.
   *
   */
  void basic_sector_decomposition::strategy_B(const integration_data & global_data, 
					      std::vector<integrand> & lst_integrand_in, std::vector<integrand> & lst_integrand_out)
  {

         std::vector<std::string> lines;
                for (int i = 0; i < lst_integrand_in.size(); i++) {
                     std::ostringstream s;
                     s << "\n(" << getIndex(global_data.list_feynman_parameter[i]);
                     s << ",(";
                     for (int j = 0; j < lst_integrand_in.size(); j++) {
                         if (i != j) {
                             s << getIndex(global_data.list_feynman_parameter[j]);
                              if (!(j == lst_integrand_in.size() - 2 && i == lst_integrand_in.size() - 1 )) {
                              if (!(j == lst_integrand_in.size() - 1)) {
                                                       s << ",";
                                                       }
                                                   }
                         }

                     }
                     s << ",),)";
                     lines.push_back(s.str());
                }

       std::cout << "sector_strings.append(\"\"\"";

    while ( lst_integrand_in.size() > 0 )
      {
	bool flag_continue = true;

	// get the last integral
	integrand local_integrand = lst_integrand_in.back();
	lst_integrand_in.pop_back();

	std::string line = lines.back();
	lines.pop_back();

	size_t n = local_integrand.index.size();
	size_t r = local_integrand.polynomial_list.size();

	// ----------------------------------------------------
	// one variable factorizes
	multi_index_ordered index_k(n,1);
	for( index_k.init(); !index_k.overflow(); index_k++)
	  {
	    for (int j=0; j<r; j++) 
	      {
		if ( local_integrand.check_zero(j, index_k, global_data) ) 
		  {
		    // can factor the k'th variable from the j'the polynomial
		    size_t k = index_k[0];
		    local_integrand.polynomial_list[j] = normal(local_integrand.polynomial_list[j]
								/global_data.list_feynman_parameter[local_integrand.index[k]]);

		    (local_integrand.weight[k]).sum_up( local_integrand.poly_exponent[j] );

		    flag_continue = false;

		  } // check_zero
	      } // j=0; j<r;
	  } // loop over index_k

	// some variables were factorized
	if ( !flag_continue) {
		    lines.push_back(line);

	lst_integrand_in.push_back(local_integrand);
		}

	// ----------------------------------------------------
	// decompose subsector, all polynomials at the same time
	if ( flag_continue )
	  {
	    std::vector<size_t> index_set;

	    if ( local_integrand.get_S_strategy_B(index_set, global_data) )
	      {
		for (int k=0; k<index_set.size(); k++)
		  {
		  		    integrand sub_integrand = local_integrand;

		    std::ostringstream os;
		    os << line << ", (" << getIndex(global_data.list_feynman_parameter[sub_integrand.index[index_set[k]]]) << ",(";
		    for (int kk=0; kk < index_set.size(); kk++) {
		      if (kk != k) {
            os << getIndex(global_data.list_feynman_parameter[sub_integrand.index[index_set[kk]]]) << ",";
		      }
		    }
		    os << "),)";
		    lines.push_back(os.str());


		    sub_integrand.create_subsector(k, index_set, global_data);

		    lst_integrand_in.push_back(sub_integrand);
		  } // k=0; k<index_set.size();

		flag_continue = false;
	      } // flag_create_subsectors
	  }

	// ----------------------------------------------------
	// no zeros

	if ( flag_continue ) {
     std::cout << line;
	lst_integrand_out.push_back(local_integrand);
	}
      }

   std::cout << "\n\"\"\")";
  }

  /**
   *
   * Iterated sector decomposition: Strategy C.
   *
   */
  void basic_sector_decomposition::strategy_C(const integration_data & global_data, 
					      std::vector<integrand> & lst_integrand_in, std::vector<integrand> & lst_integrand_out)
  {
    bool debug = false;


    std::vector<std::string> lines;
                for (int i = 0; i < lst_integrand_in.size(); i++) {
                     std::ostringstream s;
                     s << "\n(" << getIndex(global_data.list_feynman_parameter[i]);
                     s << ",(";
                     for (int j = 0; j < lst_integrand_in.size(); j++) {
                         if (i != j) {
                             s << getIndex(global_data.list_feynman_parameter[j]);
                              if (!(j == lst_integrand_in.size() - 2 && i == lst_integrand_in.size() - 1 )) {
                              if (!(j == lst_integrand_in.size() - 1)) {
                                                       s << ",";
                                                       }
                                                   }
                         }

                     }
                     s << ",),)";
                     lines.push_back(s.str());
                }

       std::cout << "sector_strings.append(\"\"\"";

    while ( lst_integrand_in.size() > 0 ) 
      {
	if (debug) std::cout << "stack " << lst_integrand_in.size() << std::endl;

	bool flag_continue = true;

	// get the last integral
	integrand local_integrand = lst_integrand_in.back();
	lst_integrand_in.pop_back();

	std::string line = lines.back();
	lines.pop_back();

	size_t n = local_integrand.index.size();
	size_t r = local_integrand.polynomial_list.size();

	if (debug)
	  {
	    for (int j=0; j<r; j++) std::cout << local_integrand.polynomial_list[j] << std::endl;
	  }

	// ----------------------------------------------------
	// one variable factorizes
	multi_index_ordered index_k(n,1);
	for( index_k.init(); !index_k.overflow(); index_k++)
	  {
	    for (int j=0; j<r; j++) 
	      {
		if ( local_integrand.check_zero(j, index_k, global_data) ) 
		  {
		    // can factor the k'th variable from the j'the polynomial
		    size_t k = index_k[0];
		    local_integrand.polynomial_list[j] = normal(local_integrand.polynomial_list[j]
								/global_data.list_feynman_parameter[local_integrand.index[k]]);

		    (local_integrand.weight[k]).sum_up( local_integrand.poly_exponent[j] );

		    flag_continue = false;

		    if (debug) std::cout << "factor " << global_data.list_feynman_parameter[local_integrand.index[k]] << " from polynom " << j << std::endl;

		  } // check_zero
	      } // j=0; j<r;
	  } // loop over index_k

	// some variables were factorized
	if ( !flag_continue) {
		 lines.push_back(line);

	lst_integrand_in.push_back(local_integrand);}

	// ----------------------------------------------------
	// decompose subsector, all polynomials at the same time
	if ( flag_continue )
	  {
	    std::vector<size_t> index_set;

	    if ( local_integrand.get_S_strategy_C(index_set, global_data) )
	      {
		if ( debug )
		  {
		    std::cout << "S: " ;
		    for (int k=0; k<index_set.size(); k++) std::cout << global_data.list_feynman_parameter[local_integrand.index[index_set[k]]] << " ";
		    std::cout << std::endl;
		  }

		for (int k=0; k<index_set.size(); k++)
		  {
		    if (debug) std::cout << global_data.list_feynman_parameter[local_integrand.index[index_set[k]]] << " ";
		    integrand sub_integrand = local_integrand;

  std::ostringstream os;
		    os << line << ", (" << getIndex(global_data.list_feynman_parameter[sub_integrand.index[index_set[k]]]) << ",(";
		    for (int kk=0; kk < index_set.size(); kk++) {
		      if (kk != k) {
            os << getIndex(global_data.list_feynman_parameter[sub_integrand.index[index_set[kk]]]) << ",";
		      }
		    }
		    os << "),)";
		    lines.push_back(os.str());

		    sub_integrand.create_subsector(k, index_set, global_data);

		    lst_integrand_in.push_back(sub_integrand);
		  } // k=0; k<index_set.size();

		flag_continue = false;
		if (debug) std::cout << std::endl;
	      } // flag_create_subsectors
	  }

	// ----------------------------------------------------
	// no zeros

	if ( flag_continue ) {
	     std::cout << line;
	lst_integrand_out.push_back(local_integrand);
	}

	if ( debug && flag_continue ) std::cout << "done." << std::endl;

	if (debug) std::cout << std::endl;
      }

   std::cout << "\n\"\"\")";
  }

  /**
   *
   * Iterated sector decomposition: Strategy X, one polynomial at the time.
   *
   */
  void basic_sector_decomposition::strategy_X(const integration_data & global_data, 
					      std::vector<integrand> & lst_integrand_in, std::vector<integrand> & lst_integrand_out)
  {

    std::vector<std::string> lines;
                for (int i = 0; i < lst_integrand_in.size(); i++) {
                     std::ostringstream s;
                     s << "(" << getIndex(global_data.list_feynman_parameter[i]);
                     s << ",(";
                     for (int j = 0; j < lst_integrand_in.size(); j++) {
                         if (i != j) {
                             s << getIndex(global_data.list_feynman_parameter[j]);
                              if (!(j == lst_integrand_in.size() - 2 && i == lst_integrand_in.size() - 1 )) {
                              if (!(j == lst_integrand_in.size() - 1)) {
                                                       s << ",";
                                                       }
                                                   }
                         }

                     }
                     s << ",),)";
                     lines.push_back(s.str());
                }

       std::cout << "sector_strings.append(\"\"\"";
        bool first = true;


    while ( lst_integrand_in.size() > 0 ) 
      {
	bool flag_continue = true;

	// get the last integral
	integrand local_integrand = lst_integrand_in.back();
	lst_integrand_in.pop_back();

	std::string line = lines.back();
    	lines.pop_back();

	size_t n = local_integrand.index.size();
	size_t r = local_integrand.polynomial_list.size();

	// ----------------------------------------------------
	// one variable factorizes
	multi_index_ordered index_k(n,1);
	for( index_k.init(); !index_k.overflow(); index_k++)
	  {
	    for (int j=0; j<r; j++) 
	      {
		if ( local_integrand.check_zero(j, index_k, global_data) ) 
		  {
		    // can factor the k'th variable from the j'the polynomial
		    size_t k = index_k[0];
		    local_integrand.polynomial_list[j] = normal(local_integrand.polynomial_list[j]
								/global_data.list_feynman_parameter[local_integrand.index[k]]);

		    (local_integrand.weight[k]).sum_up( local_integrand.poly_exponent[j] );

		    flag_continue = false;

		  } // check_zero
	      } // j=0; j<r;
	  } // loop over index_k

	// some variables were factorized
	if ( !flag_continue) {
   		lines.push_back(line);

	lst_integrand_in.push_back(local_integrand);
		}

	// ----------------------------------------------------
	// decompose subsector, one polynomial at the time
	if ( flag_continue )
	  {
	    for (int j=0; j<r; j++) 
	      {
		for (int m=2; m<=n; m++)
		  {
		    multi_index_ordered index_set(n,m);
		    for( index_set.init(); !index_set.overflow(); index_set++)
		      {
			if ( flag_continue )
			  {
			    if ( local_integrand.check_zero(j, index_set, global_data) )
			      {
				// create m new subsectors
				for (int k=0; k<m; k++)
				  {
				    integrand sub_integrand = local_integrand;

		    std::ostringstream os;
		    os << line << ", (" << getIndex(global_data.list_feynman_parameter[sub_integrand.index[index_set[k]]]) << ",(";
		    for (int kk=0; kk < index_set.size(); kk++) {
		      if (kk != k) {
            os << getIndex(global_data.list_feynman_parameter[sub_integrand.index[index_set[kk]]]) << ",";
		      }
		    }
		    os << "),)";
		    lines.push_back(os.str());
                    				    sub_integrand.create_subsector(k, index_set, global_data);

				    lst_integrand_in.push_back(sub_integrand);
				  } // k=0; k<m;

				flag_continue = false;
			      } // check_zero
			  } // flag_continue
		      } // loop over index_set
		  } // m=2; m<=n;
	      } // j=0; j<r;
	  }
	
	// ----------------------------------------------------
	// no zeros

	if ( flag_continue ) {
            		    std::cout << line << "\n";

	    lst_integrand_out.push_back(local_integrand);
	    if (first) {
	        first = false;
	    }
	}
      }

    std::cout << "\"\"\")";
  }

  /**
   *
   * Iterated sector decomposition: Strategy Y, all polynomials at the same time.
   *
   */
  void basic_sector_decomposition::strategy_Y(const integration_data & global_data, 
					      std::vector<integrand> & lst_integrand_in, std::vector<integrand> & lst_integrand_out)
  {
    while ( lst_integrand_in.size() > 0 ) 
      {
	bool flag_continue = true;

	// get the last integral
	integrand local_integrand = lst_integrand_in.back();
	lst_integrand_in.pop_back();

	size_t n = local_integrand.index.size();
	size_t r = local_integrand.polynomial_list.size();

	// ----------------------------------------------------
	// one variable factorizes
	multi_index_ordered index_k(n,1);
	for( index_k.init(); !index_k.overflow(); index_k++)
	  {
	    for (int j=0; j<r; j++) 
	      {
		if ( local_integrand.check_zero(j, index_k, global_data) ) 
		  {
		    // can factor the k'th variable from the j'the polynomial
		    size_t k = index_k[0];
		    local_integrand.polynomial_list[j] = normal(local_integrand.polynomial_list[j]
								/global_data.list_feynman_parameter[local_integrand.index[k]]);

		    (local_integrand.weight[k]).sum_up( local_integrand.poly_exponent[j] );

		    flag_continue = false;

		  } // check_zero
	      } // j=0; j<r;
	  } // loop over index_k

	// some variables were factorized
	if ( !flag_continue) lst_integrand_in.push_back(local_integrand);

	// ----------------------------------------------------
	// decompose subsector, all polynomials at the same time
	if ( flag_continue )
	  {
	    bool flag_create_subsectors = false;
	    int fact_degree = 0;
	    int best_m;
	    multi_index_ordered best_index_set;

	    for (int m=2; m<=n; m++)
	      {
		multi_index_ordered index_set(n,m);

		for( index_set.init(); !index_set.overflow(); index_set++)
		  {
		    if ( local_integrand.check_zero(index_set, global_data) )
		      {
			int actual_fact_degree = local_integrand.get_ldegree(index_set, global_data);

			if ( actual_fact_degree > fact_degree )
			  {
			    fact_degree = actual_fact_degree;
			    best_m = m;
			    best_index_set = index_set;
			    flag_create_subsectors = true;
			  } // actual_fact_degree > fact_degree

		      } // check_zero
		  } // loop over index_set
	      } // m=2; m<=n;

	    if ( flag_create_subsectors )
	      {
		// create m new subsectors
		for (int k=0; k<best_m; k++)
		  {
		    integrand sub_integrand = local_integrand;

		    sub_integrand.create_subsector(k, best_index_set, global_data);

		    lst_integrand_in.push_back(sub_integrand);
		  } // k=0; k<best_m;

		flag_continue = false;
	      } // flag_create_subsectors
	  }

	// ----------------------------------------------------
	// no zeros

	if ( flag_continue ) lst_integrand_out.push_back(local_integrand);
      }
  }

  /**
   *
   * Generate the subtraction terms.
   *
   */
  void basic_sector_decomposition::generate_subtraction_terms(const integration_data & global_data, 
							      std::vector<integrand> & lst_integrand_in, std::vector<GiNaC::ex> & lst_integrand_out)
  {
    while ( lst_integrand_in.size() > 0 ) 
      {
	// get the last integral
	integrand local_integrand = lst_integrand_in.back();
	lst_integrand_in.pop_back();

	GiNaC::ex weight_factor = local_integrand.get_weight_factor(global_data);
	GiNaC::ex polynomials   = local_integrand.get_polynomials(global_data);
	GiNaC::ex extra_factor  = local_integrand.get_extra_factor();

	GiNaC::ex poly_and_extra = extra_factor*polynomials; 
	GiNaC::ex res            = poly_and_extra;

	// index_set indicates in which variables subtractions are needed,
	// subtr_vec indicates to which order,
	// m is the size of index_set
	std::vector<size_t> subtr_vec;
	std::vector<size_t> index_set;

	int m = local_integrand.get_subtraction_vector(subtr_vec,index_set);

	// loop over all subsets of index_set
	for (int j=1; j<=m; j++)
	  {
	    multi_index_ordered subtraction_set = multi_index_ordered(m,j);
	    for ( subtraction_set.init(); !subtraction_set.overflow(); subtraction_set++)
	      {
		// local_index_set is a subset of index_set
		// no subtractions in variables, which are not in local_index_set,
		// this corresponds to no derivative in this variable.
		std::vector<size_t> local_subtr_vec;
		std::vector<size_t> local_index_set;
		for (int k=0; k<j; k++)
		  {
		    local_subtr_vec.push_back(subtr_vec[subtraction_set[k]]);
		    local_index_set.push_back(index_set[subtraction_set[k]]);
		  }

		// subtractions
		multi_index_counter_indv counter = multi_index_counter_indv(local_subtr_vec, j);
		for( counter.init(); !counter.overflow(); counter++)
		  {
		    // subtract
		    res += local_integrand.taylor_term(poly_and_extra, counter, local_index_set, global_data);

		    // and add again
		    integrand new_integrand = integrand(local_integrand.index,local_integrand.weight,poly_and_extra);
		    new_integrand.integrated_taylor_term(counter, local_index_set, global_data);
		    lst_integrand_in.push_back( new_integrand );
		  }

	      } // loop subtraction_set
	  } // j=1; j<=m;

	lst_integrand_out.push_back( weight_factor*res );
      }
  }

    /**
     *
     * Expand in epsilon.
     *
     */
    void basic_sector_decomposition::expand_in_epsilon(const integration_data & global_data,
  						     std::vector<GiNaC::ex> & lst_integrand_in, GiNaC::ex & res_out)
    {
      res_out = 0;


      GiNaC::symbol eps = GiNaC::ex_to<GiNaC::symbol>(global_data.epsilon);
      int order = global_data.order;

      for (int j=0; j<lst_integrand_in.size(); j++)
        {
  	GiNaC::ex temp = series_to_poly(lst_integrand_in[j].series(eps,order+1));
  	temp = temp.expand().collect(eps).coeff(eps,order);
  	res_out += temp;
        }
    }

  /**
   *
   * Expand in epsilon.
   *
   */
  void basic_sector_decomposition::expand_in_epsilon(const integration_data & global_data, 
						     std::vector<GiNaC::ex> & lst_integrand_in, GiNaC::ex & res_out, std::vector<GiNaC::ex> & res_out2)
  {
    res_out = 0;


    GiNaC::symbol eps = GiNaC::ex_to<GiNaC::symbol>(global_data.epsilon);
    int order = global_data.order;

    for (int j=0; j<lst_integrand_in.size(); j++)
      {
	GiNaC::ex temp = series_to_poly(lst_integrand_in[j].series(eps,order+1));
	temp = temp.expand().collect(eps).coeff(eps,order);
	res_out2.push_back(temp);
	res_out += temp;
      }
  }

  /**
   *
   * Numerical integration.
   *
   */
  void basic_sector_decomposition::numerical_integration(const integration_data & global_data, const monte_carlo_parameters & mc_parameters, GiNaC::ex res_in, monte_carlo_result & res_out, int verbose_level)
  {
    res_out = intnum( res_in, global_data.list_feynman_parameter, mc_parameters, verbose_level );
  }

} // namespace sector_decomposition

